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The  major  results  of  this  research  project  are  as  follows;  (1)  A  scheme  was  devised 
to  model  the  physical  phenomena  of  plastic  deformation,  heat  transfer,  and  phase 
transformation  in  a  unified  way.  (2)  This  scheme  was  incorporated  into  a  finite  element 
method  so  that  physical  problems  involving  the  coupling  of  all  three  phenomena  could  be 
analyzed.  (3)  Calculations  were  performed  for  two  prototype  problems:  quenching  of  a 
ball  bearing  and  of  a  cylindrical  Jominy  specimen.  (4)  Two  papers  were  written.  The  first 
is  titled  "Selective  Reduced  Integration  for  the  Lagrangian  Formulation,"  and  it  has  been 
sent  to  ARO  previously.  The  second  has  the  same  title  as  this  research  project  and  is 
included  here  as  the  main  component  of  the  final  report. 


Dr.  Daniel  C.  Peirce  was  the  principal  investigator  for  this  project.  From  time  to 


time,  he  consulted  with  Dr.  Peter  D.  Hilton  and  Dr.  John  L.  O’Brien. 
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1.  Introdpction. 

Numerical  modeling  of  materials  processing  has  expanded  considerably  in  recent 
years.  Greater  computing  power  and  broader  understanding  of  constitutive  behavior  have 
made  finite  element  calculations  of  complex  deformations  both  more  practical  and  more 
accurate.  In  addition,  heat  transfer  and  other  thermal  effects  are  now  routinely  included 
in  analyses  of  such  processes  as  rolling  and  forging  [1—3].  Models  of  quenching  and  heat 
treatment,  however,  have  developed  separately  from  efforts  on  processes  involving 
materials  forming  [4—5],  no  doubt  due  to  the  relatively  lesser  roles  of  heat  transfer  and 
phase  transformation  in  the  latter,  together  with  the  relatively  greater  extent  of  plastic 
flow.  The  aim  of  this  paper  is  to  establish  a  modeling  framework  wherein  the  effects  of 
heat  transfer  and  phase  change  may  be  coupled  with  plastic  deformation  in  a  unified  way. 
While  specific  attention  here  will  be  focused  on  examples  of  quenching  procedures,  the 
framework  developed  should  broaden  the  capabilities  of  modeling  for  a  wide  range  of 
materials  processes. 

A  brief  review  of  previous  work  which  has  proceeded  in  this  direction  indicates  the 
need  for  the  type  of  unifying  approach  that  will  be  proposed  here.  Inoue  and  Raniecki  [6] 
presented  a  table  of  research  papers  which  had  dealt  with  determination  of  quenching 
stresses,  but  only  a  few  of  the  studies  referenced  in  their  table  combined  all  three  effects  of 
phase  change,  thermal  stress,  and  plasticity.  Inoue  and  Raniecki  built  on  these  works  to 
develop  a  careful  treatment  of  the  quenching  problem,  with  special  attention  given  to  the 
plasticity  formulation.  Nevertheless,  their  analyses  of  the  heat  transfer  and  corresponding 
phase  transformation  history  is  completed  prior  to  and  independently  of  their  solution  of 
the  plasticity  problem  (which  predicts  residual  stresses). 

Subsequent  work  by  Sjostrom  [7]  and  Hildenwall  and  Ericsson  [8]  also  included  the 
effects  of  plastic  flow,  but  the  heat  transfer  and  phase  transformation  calculation  remained 
decoupled  fiom  the  stress  analysis.  The  studies  by  Fischer  et  al  [9,10]  kept  the  stress 
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analysis  separate  too,  but  they  took  special  care  to  account  for  the  temperature 
dependence  of  material  properties  for  the  heat  transfer  and  phase  change  calculation.  None 
of  these  methodologies,  though,  includes  the  effect  of  latent  heat  of  transformation, 
heating  due  to  plastic  work,  and  the  stress  influence  on  phase  changes.  While  these  effects 
may  be  small  in  many  situations,  a  flexible  and  unified  approach  should  allow  for  their 
presence. 

Inoue  and  Wang  [11]  recognized  the  completeness  of  this  coupling  which  exists 
among  heat  transfer,  phase  change,  and  plastic  flow.  Denis  et  al  [4]  also  took  account  of 
the  stress  effects  on  martensitic  transformation,  and  emphasized  their  treatment  of  the 
transformation  distortion  as  additional  plastic  strain.  These  papers  provide  careful 
analysis  of  the  phase  change  and  thermal  effects,  but  these  are  coupled  with  the  plastic 
flow  by  using  a  simple  incremental  approach  in  which  the  heat  transfer  and  stiffness 
equations  are  solved  in  alternating  fashion. 

This  incremental  approach  to  coupled  analysis  has  been  used  often  for  research  on 
thermomechanical  problems  which  do  not  involve  phase  transformations.  For  example, 
Dawson  and  his  coworkers  have  developed  finite  element  programs  to  model 
thermomechanical  response  in  such  applications  as  salt  mines  and  forming  of  porous 
materials  [12—14].  Chandra  and  Mukherjee  [2]  developed  a  coupled  formulation  for  forming 
processes  with  heat  transfer,  and  they  used  it  to  solve  an  extrusion  problem.  In  contrast 
to  these  approaches,  Rebelo  and  Kobayashi  [1]  developed  a  procedure  for  coupled  analysis 
which  employs  an  iterative  scheme.  The  iterations  serve  to  maintain  close  contact  between 
the  thermal  and  mechanical  solutions.  Absent  from  their  constitutive  equations,  however, 
are  elastic  deformations;  these  are  important  in  many  situations,  such  as  when  an 
accurate  determination  of  residual  stresses  is  required. 

The  modeling  framework  presented  here  will  attempt  to  bring  together  the  best 
features  of  both  the  quenching  and  the  metal  forming  research  represented  in  the  foregoing. 
The  intent  is  not  to  give  the  final  answer  on  which  constitutive  models  and  which  phase 
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transformation  laws  are  the  correct  ones  for  solving  materials  processing  problems.  Rather, 
the  goal  here  is  to  formulate  equations  which  are  general  enough  to  accommodate  readily 
any  or  all  of  the  key  phenomena  which  occur  in  thermomechanical  problems  with  phase 
changes.  The  pertinent  aspects  of  these  phenomena  will  be  further  discussed  below,  but 
reference  is  made  here  to  the  chart  presented  by  Inoue  and  Wang  [11]  and  reproduced  with 
slight  variation  in  Figure  1.  This  chart  shows  that  temperature  and  stress  affect  each 
other,  that  stress  and  phase  change  affect  each  other,  and  that  phase  change  and 
temperature  affect  each  other.  The  equations  to  be  developed  here  will  be  able  to  model  all 
these  interactions,  and  they  will  provide  for  great  flexibility  in  the  descriptions  of 
plasticity  and  phase  changes.  The  formulation  established  below  has  been  applied  to 
example  problems  of  the  quenching  of  a  ball  bearing  and  a  Jominy  end  quench.  The  details 
of  these  and  future  applications  are  reserved  for  future  reports. 


2.  Modeling  Background. 

This  section  discusses  large  deformation  mechanics,  phase  transformation,  and 
heat  transfer  in  turn.  The  separate  models  for  these  phenomena  will  provide  the  basis  for 
the  unified  coupling  to  be  described  in  Section  3. 

2.1  Large  deformations  and  plasticity.  The  ultimate  use  for  the  methodology 
proposed  here  is  its  incorporation  into  a  finite  element  framework.  In  such  a  framework, 
incremental  equilibrium  is  enforced  by  relying  on  the  incremental  principle  of  virtual  work 


I  t  •  dV 
Vo 


I  T  .  tfudS. 
So 


(2.1) 


Here,  the  nominal  stress  t  and  deformation  gradient  F  are  calculated  with  respect  to  the 
reference  configuration,  which  has  volume  Vq  and  surface  Sq.  In  many  cases,  the 
reference  configuration  corresponds  to  the  initial  undeformed  state  of  the  body.  The  vector 
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T  is  the  nominal  traction  vector,  and  u  is  the  displacement  vector.  Constitutive 
relations  are  often  expressed  conveniently  in  terms  of  the  Kirchhoff  stress  r,  which  is 
related  to  the  nominal  stress  t  and  Cauchy  stress  c  by 

r  =  F  •  t  =  Ja,  (2.2) 

where  J,  the  Jacobian  of  the  deformation,  is  the  determinant  of  F. 

For  material  models  which  incorporate  dastic  behavior,  incremental  response  is 
based  on  the  relationship 

?■  =  L  :  De,  (2.3) 

where  r  is  the  Jaumann  rate  of  Kirchhoff  stress,  L  is  the  tensor  of  dastic  moduli,  and 
De  is  the  elastic  portion  of  the  deformation  rate  tensor  D,  which  is  the  symmetric  part  of 
F  •  F‘^  Further  details  of  the  Jaumann  rate  and  the  kinematics  behind  equation  (2.3) 
may  be  had  from  a  variety  of  sources,  including  Needleman  [15,16]  and  Anand  ti  al  [17]. 
For  the  kinds  of  problems  under  consideration  here,  the  deformation  rate  D  is  made  up  of 
several  distinct  constituents  in  addition  to  D«: 

D  =  D«  +  Dvp  +  Dtr  +  Dte,  (2.4) 

where  D^p  represents  the  contribution  from  viscoplasticity,  contains  the  effects  of 
phase  transformation,  and  incorporates  the  influence  of  thermal  expansion  and 
contraction.  Each  of  these  terms  must  be  expressed  in  such  a  way  that  (2.3)  may  be 
inserted  into  (2.1)  and  stiffness  equations  derived.  This  subsection  focuses  on  D'^p,  while 
and  will  be  discussed  in  connection  vnth  phase  transformation  and  heat  transfer, 
respectively. 
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A  viscoplastic  formulation  is  used  here  because  strain— rate  dependence  and  creep 
play  important  roles  in  many  materials  processes.  With  some  loss  of  generality,  attention 
will  be  restricted  to  materials  which  harden  isotropically,  and  for  which  the  plastic 
deformation  rate  D^p  is  always  aligned  with  the  Kirchhoff  stress  deviator  s.  For 
rate— independent  materials,  these  assumptions  are  consistent  with  J2  flow  theory. 

Potential  variations  from  these  restrictions  include  pressure-sensitive  plastic  flow  and 
kinematic  hardening  behavior.  (Indeed,  the  former  will  play  a  role  here  due  to  equation 
(2.10)  introduced  later.)  There  is  no  reason  that  such  behavior  cannot  be  included  in  what 
follows,  and  reference  is  made  to  Anand  et  al  [17]  and  Peirce  et  al  [18],  where  further 
details  on  modeling  these  and  other  effects  may  be  found.  The  presentation  here,  however, 
becomes  involved  enough  even  with  the  various  assumptions  put  forth  above. 

The  Kirchhoff  stress  deviator  and  related  quantities  are  defined  as  follows: 

s  s  tI  $  |«;a,  (2.5a) 

p  a  P  a  L ;  p.  (2.5b) 

Here  re  represents  the  effective  Kirchhoff  stress,  and  p  and  P  are  tensors  which  prove 
convenient  in  the  subsequent  development.  It  follows  from  (2.5b)  that  p  :  p  =  3/2.  The 
viscoplastic  deformation  rate  may  now  be  written 

Dvp  =  £vpp^  (2.6) 

where  the  accumulated  viscoplastic  strain  is  defined  as 
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A  tremendous  research  effort  has  been  expended  on  determination  of  useful  functional 
forms  for  c'^p,  which  may  in  general  depend  on  stress,  temperature,  plastic  work,  and 
several  state  variables.  Here  again,  attention  wiU  be  restricted  to  a  limited  subset  of  the 
general  case,  to  be  represented  as  follows: 

cvp  =  f(re,  fvp,  T,  Pi),  (2.8) 

where  T  is  temperature,  and  each  pi  here  represents  the  volume  fraction  of  phase  i  of 
the  material.  The  development  which  follows  could  be  carried  out  equally  well  with 
different  choices  for  (2.8);  examples  include  the  models  proposed  by  Merzer  and 
Bodner  [19]  and  Anand  [20]. 

At  this  point,  the  form  (2.6)  with  (2.8)  specified  is  enough  to  provide  the  plasticity 
input  to  the  incremental  virtual  work  principle  (2.1).  The  functional  expressions  typically 
associated  with  (2.8)  give  rise  to  stiff  response  when  time  integration  is  attempted. 
Moderation  of  this  stiff  behavior  is  a  key  reason  for  the  development  of  the  integration 
techniques  to  be  developed  in  Section  3. 

2.2  Phase  transformations.  The  distortion  induced  by  a  phase  transformation, 
represented  here  by  arises  from  volume  changes  associated  with  the  growth  of  a  new 
phase.  The  volume  increase  which  typically  accompanies  a  martensitic  transformation  is 
an  excellent  example  of  this  phenomenon.  The  magnitudes  of  such  volume  changes  are 
generally  well-characterized,  and  will  be  denoted  here  by  3bi,  on  a  percentage  basis,  for 
appearance  of  each  phase  i.  Therefore, 


Dtr  =  I  pi  Pi-  (2.9) 

The  analogy  between  bi  and  the  coefficient  of  thermal  expansion  is  worth  noting,  to  be 
completely  general,  it  would  be  better  to  write  (2.9)  with  two  subscripts  on  each  character 
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to  correspond  to  a  specific  transformation  (from  phase  i  to  phase  j,  say).  It  is  often  the 
case,  however,  that  a  specific  phase  is  formed  in  only  one  way;  in  addition,  the  number  of 
phases  worth  considering  in  many  applications  will  likely  be  three  or  fewer,  so  that  (2.9) 
should  not  be  the  source  of  any  confusion. 

While  accurate  values  for  bi  are  generally  quite  accessible,  the  nucleation  and 
growth  of  a  new  phase  can  proceed  in  any  number  of  ways,  and  suitably  realistic  models 
for  the  evolution  of  the  phases  pi  are  complex.  For  the  purposes  of  this  work,  it  will  be 
assumed  that  the  pi  may  be  represented  in  functional  form  as 


Pi  =  Pi  (U,  Th,  T,  t), 


(2.10) 


where  T  is  temperature  and  t  is  time;  of  course,  the  pi  must  sum  to  unity.  The 
effective  stress  re  has  been  defined  in  (2.5a),  while  rh  is  roughly  thrice  the  negative 
pressure  — p,  that  is 


n  =  I:  r  =  J  (I :  a)  =  -3pJ. 


(2.11) 


Although  situations  no  doubt  exist  where  (2.10)  is  too  restrictive,  such  as  when  chemistry 
or  grain  size  come  into  play,  these  factors  are  usually  fixed  during  a  given  process. 
Equation  (2.10)  does  account  for  some  of  the  key  driving  forces  for  phase  changes:  stress, 
temperature,  and  time. 

The  transformation  of  austenite  to  pearlite  and  martensite  during  quenching  of  steel 
furnishes  a  useful  example  of  possible  choices  for  specification  of  (2.10).  These  phases  will 
be  denoted  here  by  pa,  Pp,  and  pn  respectively,  with  Pa  =  1  —  Pp  —  Pm.  For  pearlitic 
growth,  an  equation  along  the  lines  of  the  sigmoidal  curve  proposed  by  Johnson  and 
Mehl  [21]  may  be  used: 
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Pp  =  1 -exp[-a(T)t“]. 


(2.12) 


The  form  of  a(T)  depends  on  the  nudeation  and  growth  rates  in  the  material  at  hand,  but 
its  explidt  variation  with  temperature  T  allows  it  to  account  for  the  knee  in  the 
time— temperature-transformation  curve.  The  exponent  n  is  related  to  the  shape  of  the 
growing  phase.  More  detailed  possibilities  for  (2.12)  are  discussed  at  length  by  Inoue  and 
Raniecki  [6].  For  the  martensite  transformation,  an  expression  of  the  form 

Pb  =  (1 -Pp)|l-exp[^(T-Ms)]|, 

valid  only  when  T  is  less  than  the  martensite  start  temperature  Ms,  is  discussed  by 
Inoue  and  Raniecki  [6]  and  Hakberg  and  Hogberg  [5],  among  many  others.  The  parameter 
0  is  presumed  constant  here. 

Equations  (2.12)  and  (2.13)  show  no  dependence  on  stress,  but  this  may  be 
introduced  easily  by  allovdng  a(T)  or  Ms  to  change  with  pressure  or  effective  stress. 
Because  formation  of  martensite  results  in  an  instantaneous  volume  increase,  the  Mg 
temperature  is  effectively  increased  under  hydrostatic  tension,  and  decreased  under 
hydrostatic  pressure.  This  kind  of  effect  may  be  easily  incorporated  in  (2.13).  While 
particular  cases  may  require  greater  detail  than  is  contained  in  the  above  expressions  for 
pearlite  and  martensite  growth,  the  latter  nevertheless  indicate  how  key  aspects  of  these 
phase  transformations  may  be  modeled. 

2.3  Heat  transfer.  The  first  role  of  heat  transfer  is  to  drive  thermal  stresses  via  D<^e 
in  (2.4).  Under  the  assumption  that  each  phase  pi  has  a  coefficient  of  thermal  expansion 
ori,  a  reasonable  approximation  is 


(2.13) 


DM  =  I  (S  oi  Pi)  T. 
i 


(2.14) 
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If  appropriate,  more  accurate  expressions  for  the  coefficient  of  thermal  expansion, 
parenthesized  in  (2.14),  may  be  found  in  Bobeth  and  Diener  [22]  and  Hashin  [23]. 
Temperature  also  has  influence  by  its  effect  on  viscoplasticity,  as  in  (2.8),  and  especially 
on  phase  transformation,  as  in  (2.12)  and  (2.13).  In  addition,  many  material  properties, 
including  elastic  moduli,  coefficients  of  thermal  expansion,  and  the  thermal  properties  to 
be  introduced  in  this  section,  may  vary  significantly  with  temperature. 

For  all  these  reasons,  the  heat  transfer  equation  must  be  solved,  together  with  the 
equilibrium  equation  (2.1),  to  provide  a  complete  description  of  many  materials  processes. 
The  form  of  the  heat  transfer  equation  employed  here  is 

/iCpT  =  !.(«VT)  +  q,  (2.15) 

where  p  is  the  material  density,  Cp  is  the  specific  heat,  k  is  thermal  conductivity,  and 
q  represents  all  heat  source  terms.  For  multiphase  materials  in  which  the  phases  have 
different  thermal  properties,  weighted  averages  of  these  properties  may  be  substituted  in 
(2.15)  if  the  property  variations  among  the  phases  are  not  too  great.  Because  equation 
(2.15)  will  be  used  to  model  heat  flow  in  a  deforming  material,  the  meaning  of  the  gradient 
must  be  considered  carefully.  As  discussed  by  Lemonds  and  Needleman  [24],  heat 
conduction  depends  at  root  on  the  spacing  between  atoms,  which  is  not  affected  by 
plasticity,  but  only  by  elastic  straining.  Therefore,  if  the  gradient  in  (2.15)  is  taken  with 
respect  to  the  reference  configuration,  the  error  incurred  will  be  of  the  order  of  the  elastic 
strains,  which  are  typically  much  less  than  one  percent.  Adjustment  of  the  gradient  with 
respect  to  the  elastic  deformation  may  be  performed  if  this  additional  accuracy  is  needed, 
but  in  a  numerical  calculation,  substantial  extra  work  associated  with  splitting  the 
deformation  gradient  would  be  required.  If  significant  changes  in  atom  spacing  are 
associated  with  phase  transformation,  they  are  perhaps  best  handled  through  adjustment 
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of  the  thermal  conductivity  «i  associated  with  the  particular  phase. 

Further  influence  of  the  deforming  material  on  (2.15)  is  felt  through  the  source 
terms  q.  In  the  present  development,  these  source  terms  will  include  the  heating  due  to 
plastic  working  and  the  latent  heat  of  transformation: 

q  =  X  ^ +  S  (pb)i  pi.  (2.16) 

In  this  equation,  the  parameter  x>  which  represents  the  fraction  of  plastic  work  which 
actually  contributes  to  a  temperature  rise  in  the  body,  was  used  by  Lemonds  and 
Needleman  [24].  When  one  phase  is  transformed  into  another  phase  pi,  the  latent  heat  is 
denoted  here  by  hi,  expressed  as  energy  per  unit  mass  of  transformation  product.  With 
the  corresponding  density  pi,  (ph)i  represents  the  latent  heat  per  unit  volume  of  newly 
produced  phase  i.  The  same  remarks  about  two  subscripts  which  were  made  in  connection 
with  (2.9)  also  apply  to  (2.16),  but,  at  the  risk  of  some  confusion,  the  single  subscript  is 
used  for  brevity. 

It  is  clear  from  equations  (2.4)  and  (2.15)  and  their  ingredients  that  that  the 
thermal  and  mechanical  responses  of  many  materials  processes  are  inseparably  coupled. 

For  solution  of  these  equations,  a  numerical  approach  which  reflects  this  coupling  therefore 
has  great  appeal.  The  methodology  to  be  developed  here  attempts  to  exploit  the 
interdependence  of  the  relevant  field  quantities  in  an  incremental  scheme  which  follows  the 
true  evolution  of  these  quantities  accurately.  In  addition,  maintaining  the  closeness  of  the 
coupling  throughout  the  numerical  approach  reduces  the  tendency  to  instability  and  may 
allow  significantly  larger  time  increments  to  be  used. 

3.  Numerical  method. 

Even  with  the  various  special  assumptions  made  throughout  Section  2,  there 
remains  a  great  deal  of  complexity  in  the  models  of  the  physical  phenomena  discussed. 
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This  section  therefore  begins  with  a  summary  presentation  of  the  final  numerical  scheme, 
with  description  of  its  full  detail  and  derivation  to  follow. 

3.1  Coupling  of  the  mechanical  and  thennal  solutions.  As  noted  in  Section  2,  the 
mechanical  response  is  kept  in  equilibrium  by  the  incremental  virtual  work  principle  (2.1). 
The  method  for  implementing  (2.3)  in  (2.1)  may  be  found  elsewhere  [15,16].  The  key  task 
for  modeling  the  mechanical  response  is  thus  to  specify  (2.3)  in  detail.  As  subsequent 
manipulations  will  show,  the  Jaumann  rate  of  Kirchhoff  stress  takes  the  form 

?■  =  C  ;  D  +  Gs  P  +  Gg  I  +  (Gt  P  +  Gg  I)  T,  (3.1) 

where  Gk  =  N'^Hk,  C  is  given  by 

C  =  N-ifL-HtPP-HsPI-HglP-Hgll],  (3.2) 

and  N  and  the  Hk  for  k  =  1  to  8  are  defined  later.  At  this  point,  it  should  be  noted 

that  (3.1)  depends  on  deformation  rate  D,  temperature  rate  T,  and  the  other  terms  in 
the  middle  of  the  right-hand  side,  which  may  be  specified  completely  in  terms  of  the 
current  state.  These  terms  are  analogous  to  the  residual  loading  terms  which  appear  in 
many  plasticity  formulations.  In  (3.2),  unless  H2  =  Hg,  the  moduli  C  do  not  lead  to  a 
symmetric  stiffness  matrix  in  a  finite  element  formulation.  This  aberration  arises  quite 
expectedly  from  the  stress  dependence  of  the  phase  transformations  as  indicated  in  (2.10). 

It  will  be  seen  later  that  if  stresses  do  not  appear  in  (2.10),  Hg  =  Hg  =  0. 

Treatment  of  the  heat  transfer  equation  by  finite  elements  requires  that  (2.15)  be 
multiplied  by  a  variation  in  temperature  ST  and  integrated  over  a  reference  volume.  As 
in  [24],  this  results  in 
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f  p  Cp  T  dV  =  -  I  /c  (V  T)  .  (V  <5T)  dV  + 

•'Vo  •^Vo 

)  q  dV  +  )  «  (n  •  V  T)  dS. 
*'Vo  •'So 


(3.3) 


Here,  n  is  the  unit  outward  normal  to  the  body  in  the  reference  configuration.  Equation 

(3.3)  leads  to  the  following  equations  for  the  nodal  temperatures  Tj  and  their  rates  Tj 
[24]: 

Mij  Tj  +  Dij  Tj  =  Qi.  (3.4) 


Here,  the  subscripts  i  and  j  refer  to  finite  dement  nodes  and  the  matrices  represented 
by  Mij  and  Dij  are  derived  in  standard  fashion  fi-om  the  integrals  in  (3.3).  The 
following  approach  is  adopted  here  for  the  solution  of  (3.4).  For  some  0,  0  <  ^  <  1, 
equation  (3.4)  is  presumed  to  hold  at  time  t  +  0At,  where  t  is  the  time  and  At  is  the 
time  increment.  The  value  of  Tj  at  t  +  ^At  is  linearly  interpolated  from  its  values  at  t 
and  t  +  At.  As  a  result,  the  equation  may  be  multiplied  by  At  and  rearranged  to  give 

(Mij  +  ^At  Dij)  Tj(t  +  At)  = 

(Mij  -  (1  -  ^At  Dij)  Tj(t)  +  Qi(t  +  ^At),  (3.5) 

which  may  readily  be  solved  for  Tj(t  +  At).  The  nodal  values  Qi(t  +  6Ai)  arise  from 
(2.16)  and  depend  on  any  source  terms  that  may  be  present  as  well  as  on  the  current 
increment  of  mechanical  behavior  through  the  plastic  work  increment  derived  later  in 
(3.14). 

Equation  (3.4)  in  turn  is  coupled  with  the  mechanical  equilibrium  equations  for  the 
nodal  displacement  rates  Uj: 
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Kij  Uj  +  Cij  Tj  =  Fi. 


(3.6) 


The  finite  element  matrices  Kij  and  Cy  are  derived  from  inserting  (3.1)  into  (2.1).  Both 

(3.4)  and  (3.6)  depend  on  increments  of  both  the  nodal  temperatures  and  the  nodal 
displacements;  the  latter  enter  (3.4)  via  the  vector  of  source  terms  Qi  which  derive  from 
(2.16). 

The  solution  espoused  here  for  the  solution  of  the  coupled  equations  (3.4)  and  (3.6) 
is  somewhat  similar  to  that  put  forth  by  Lemonds  and  Needleman  [24].  At  a  particular 

time  t,  it  is  assumed  that  the  current  state  is  known  and  that  the  incremental  fields  Tj 

and  Uj  must  be  determined.  First,  an  estimate  of  the  temperature  increments  Tj  is 
obtained  as  follows.  At  each  nodal  point,  a  quadratic  is  fit  to  the  values  of  Tj  at  the 
current  time  t  and  the  two  previous  time  steps.  This  quadratic  is  used  to  obtain  a  value 
of  each  Tj  at  time  t  +  At,  from  which  the  temperature  increment  may  be  obtained  by 

subtracting  Tj(t).  With  these  estimates  of  the  Tj,  equation  (3.6)  may  be  solved  for  the 

Uj.  This  solution  in  turn  provides  the  information  required  to  specify  Qi(t  +  ^At)  and  to 

solve  (3.4)  for  a  presumably  more  accurate  set  of  the  Tj.  This  new  set  may  be  compared 
with  the  estimates  from  the  quadratic  fit,  and  if  the  difference  is  higher  than  a  prescribed 
tolerance,  an  iterative  procedure  may  be  undertaken,  or,  as  is  the  case  with  current 
applications  of  this  approach,  the  entire  increment  may  be  recalculated  with  a  smaller 
time  step.  To  initialize  a  calculation,  when  previous  temperature  data  are  unavailable,  it 
seems  best  to  solve  (3.4)  without  taking  account  of  the  mechanical  behavior  inputs  to  Qi; 
here  again,  the  estimates  thus  arrived  at  can  be  checked  once  (3.6)  has  been  solved  and 

(3.4)  solved  again  with  the  new  information.  The  main  drawback  of  repeated  solutions  of 

(3.4)  is  that  they  are  much  more  expensive  than  quadratic  fits.  In  quenching  processes  and 
other  situations  involving  steep  temperature  gradients,  however,  a  quadratic  fit  tends  to 
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be  less  accurate  than  direct  solution  of  (3.4)  without  the  complete  Qi,  because  the  rapid 
heat  conduction  dominates  what  is  typically  a  smaller  influence  from  the  heat  source  terms. 

Once  the  incremental  fields  Uj  and  Tj  have  been  obtained  from  (3.4)  and  (3.6), 
the  stresses,  phases,  and  other  internal  variables  may  be  updated.  Finally,  the 
displacements  Uj  and  temperatures  Tj  at  time  t  +  At  are  calculated.  At  this  point, 
t  +  At  becomes  the  current  time,  and  the  next  increment  may  begin. 

3.2  Derivation  of  the  constitutive  laws.  It  remains  to  fill  in  the  details  of  (3.1)  and 
(3.2).  The  starting  point  is  the  viscoplastic  deformation  rate  D^p  in  equations  (2.4)  and 
(2.6).  Because  typical  forms  for  (2.8)  lead  to  stiff  constitutive  equations,  the  tangent 
modulus  method  proposed  by  Peirce,  Shih,  and  Needleman  [18]  and  extended  by  Anand  et 
al  [17]  is  followed  in  the  present  development  in  order  to  alleviate  problems  associated  with 
the  stiff  equations.  The  tangent  modulus  method  begins  with  the  following  frequently 
employed  expression  for  D^p  in  the  middle  of  a  time  step: 

Dvp  =  (l-tf)Dvp(t)  +  ^D''P(t  +  At).  (3.7) 

A  simple  expansion  of  (2.6)  from  the  current  time  t  is  used  to  arrive  at  the  approximation 

Dvp  =  (f  +  f>Af)  p  +  f^At  J,  (3.8) 

where  Af,  by  virtue  of  (2.8),  is 

and  may  be  determined  from  (2.5b); 
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(3.10) 


Anand  et  al  [17]  included  a  term  corresponding  to  the  rightmost  term  in  (3.8)  to  account 
for  the  change  in  direction  of  plastic  flow  during  a  time  increment,  and  this  has  the  added 
benefit  of  improving  the  ability  of  the  scheme  to  follow  nonproportional  loading  accurately. 
Note  that  when  \  is  parallel  to  s  and  hence  to  p,  the  expression  for  ^  in  (3.10) 
vanishes. 

The  next  step  is  expression  of  (3.9)  in  such  a  way  that  stress  and  phase  increments 
are  eliminated.  From  (2.5a)  and  (2.3),  there  follows 


Are  =  P:DAt  -  PrpAfvp, 


(3.11) 


while  (2.10)  leads  to 

Api  =  I^Are  +  ^Arh  +  ll^AT  +  At.  (3.12) 

Unfortunately,  equation  (3.12)  introduces  the  additional  quantity  Ath,  and  this  in  turn 
requires  that  attention  be  shifted  to  and  in  (2.10)  and  (2.14).  Equations  (2.3) 
and  (2.4)  lead  to 


Arb  =  I :  ?•  At  = 

I :  L  :  (D  At  -  I E  bi  Api  -  I  (E  tti  pi)  AT).  (3.13) 

After  the  identification  Ae^p  =  (f  +  (?Af)  At  has  been  made,  simultaneous  solution  of 
(3.9)  and  (3.11)  through  (3.13)  is  required  to  determine  expressions  for  the  incremental 
quantities  above  in  terms  of  just  At,  DAt,  and  AT.  The  algebra  involved  in  carrying 
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out  this  solution  is  quite  cumbersome.  The  details  are  made  more  tractable  by  using  (3.12) 
only  after  it  has  been  summed  against  the  phase  change  expansion  coefGdents  bi,  and  by 
seeking  expressions  for  incremental  scalar  quantities  that  represent  sums  of  four  terms, 
each  of  which  involves  At,  P:D  At,  I:D  At,  and  AT.  The  procedure  is  not  unlike  that 
employed  to  invert  Hooke’s  law  when  using  index  notation.  The  key  results  of  these 
manipulations  may  be  summarized  as  follows: 


(f  +  0A{)  At  =  F,At  +  F2P;D  At  +  FsLD  At  +  F4AT,  (3.14) 

E  bi  Api  =  BiAt  +  BjPrD  At  +  BjIrD  At  +  B4AT.  (3.15) 

For  simplicity,  the  assumption  made  here  and  from  henceforth  is  that  the  elastic  moduli 
are  isotropic,  and  that  they  are  specified  by  a  bulk  modulus  x  and  shear  modulus  fi.  On 
this  basis,  the  expressions  for  Fi  through  F4  are  given  by 

Fi  =  /3[f+ Mt(Ai-9x7A3Ri)] 

F2  =  /3  [^At(A2  —  9«7A3R2)] 

F3  =  /?  [0At(3x7A3)] 

Fi  =  P  [5At(A4  -  9x7A3(R4  +  SaiPi))], 

where 

P  =  [1  +  0At(3/iA2-A5-27x^7A3R2)]'‘,  (3.17) 

7  =  [l  +  9xSbi|^]-S  (3.18) 


(3.16a) 

(3.16b) 

(3.16c) 

(3.16d) 


and  Ai  through  A  5  are  given  by 
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"  ?^il?h’ 

A4  =  If  + 

^5  =  (3.19e) 

Further,  the  definitions  of  Ri  through  R4  are 
R.  =  Sbi^‘,  R,  = 

These  allow  Bi  through  B4  to  be  written  as 

Bi  =  7  (Rj— S/iFiRj/,  (3.21a) 

^2  =  7  (^2“  3/4F2R2)»  (3.21b) 

B3  =  7  (S/iRs— 3^F3R2),  (3.21c) 

B4  =  7  (R4- 3/XF4R2  -  9/CR3  SftiPi).  (3.21d) 

Finally,  the  details  of  (3.1)  and  (3.2)  may  be  filled  in.  The  definition  of  N  maybe 
recognized  relatively  easily  from  (3.8)  and  (3.10): 

N  =  1  +  (3.22) 

The  definitions  of  Hi  through  Hg  are  as  follows: 
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Hi  =  NF2 

T  e 

(3.23a) 

H2  H  NF3 

(3.23b) 

H3  =  3/SB2N 

(3.23c) 

H4  H  3«B3N  - 

■»  e 

(3.23d) 

Hs  E  -NF, 

(3.23e) 

Hg  =  — 3/cBiN 

(3.23f) 

H7  =  -NF4 

(3.23g) 

Hg  5  — ^3/c  (B4  -|-  SotiPi)  N 

(3.23h) 

These  equations,  while  cumbersome,  specify  the  constitutive  law  (3.1)  completely. 

Several  assumptions  have  been  made  during  this  derivation,  but  there  are  also  several  key 
effects  which  have  not  been  written  off.  For  example,  the  important  influences  of  pressure 
on  phase  changes  has  been  included;  this  is  reflected  in  the  lack  of  symmetry  of  the  moduli 
C  in  (3.2).  The  complexity  of  equations  (3.16)  to  (3.23)  may  be  collapsed  considerably  if 
additional  assumptions  are  made  with  respect  to  equations  (2.8)  or  (2.10).  If  more 
generality  is  needed  than  these  equations  afford,  the  manipulations  summarized  above 
may  be  redone  in  similar  fashion  to  obtain  appropriate  parallels  to  (3.23). 

4.  Conclusion  and  Fatore  Work.  Perhaps  the  primary  impetus  for  development  of  the 
modeling  procedures  described  above  is  the  need  for  solutions  to  problems  which  occur 
during  quenching  processes,  such  as  adverse  residual  stresses  and  quench  cracking.  It  is 
during  such  processes  that  phase  transformations  join  thermomechanical  response  as  key 
features  of  the  material  behavior.  But  phase  changes  play  important  roles  in  other 
processes  as  well,  and  they  may  be  driven  more  by  deformation  than  by  temperature  in 
bulk  fornaing  processes  like  rolling  and  forging.  The  methodology  described  above  can 
easily  be  simplified  or  enhanced  depending  on  the  important  parameters  in  a  particular 
problem.  For  example,  the  temperature  dependence  of  elastic  moduli  may  play  a  key  role 
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in  many  processing  problems;  even  though  such  dependence  has  not  been  included  here,  it 
is  typically  smooth  and  may  be  easily  incorporated  into  the  formulation. 

The  initial  application  of  the  work  in  this  manuscript  uses  the  present  methodology 
to  analyze  two  quenching  problems;  a  ball  bearing  and  an  end  quench.  The  ball  bearing 
problem  was  analyzed  using  a  simplified  version  of  the  procedure  discussed  above.  In 
particular,  the  mechanical  effects  on  the  heat  conduction  behavior  were  neglected  so  that 
the  thermal  response  was  independent  of  the  mechanical  and  phase  responses,  which  were 
coupled.  The  end  quench  problem  did  not  make  these  simplifications,  as  it  included  all 
the  coupling  effects  described  in  the  foregoing  analysis.  The  computations  performed  for 
these  two  quenching  problems  illustrated  some  important  considerations  for  conducting 
such  numerical  calculations,  especially  those  in  which  high  gradients  are  involved,  and 
these  will  be  discussed  in  subsequent  work.  Beyond  the  solution  of  these  sample  problems, 
it  is  anticipated  that  the  present  work  will  soon  find  practical  application  to  the  solution  of 
processing  problems  in  gun  barrel  quenching. 
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FIGURE  1 


This  schematic  indicates  the  key  coupling  mechanisms  in  materials 
processing  problems. 
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